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Abstract 



We consider a two-parameter nonhermitean quantum- mechanical Hamil- 
tonian operator that is invariant under the combined effects of parity and 
time reversal transformations. Numerical investigation shows that for 
some values of the potential parameters the Hamiltonian operator sup- 
ports real eigenvalues and localized eigenfunctions. In contrast with other 
PT symmetric models which require special integration paths in the com- 
plex plane, our model is integrable along a line parallel to the real axis. 



PACS numbers: 03.65.-w, 03.65. Ge, 02.60.Lj, 11.30.Er, 12.90.+b, 34.10. +x, 
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1 Introduction 



Recent analysis of the spectra of the family of Schrodinger operators with a 
complex PT invariant potential [jjj 

has raised considerable interest on such class of operators The motivation 

for the requirement of Parity times Time Reversal (PT) (actually, complex 
conjugation) symmetry comes from a conjecture of D. Bessis relating it to the 
existence of real energy bound states (in other words, non-decaying resonances). 
From a practical point of view, PT symmetry is a way of selecting a specific 
Riemann sheet in the (— oo, 0) cut complex coordinate plane, the cut being 
necessary to cope with the x = branch point of the potential. 

The loss of hermiticity of H seems to imply complex energies for localized 
eigenstates, with Im E n ^ (i.e., levels with nonzero widths), and consequently 
certain rate of decay of any initially localized state. The loss of hermiticity may 
be mediated, e.g., by a spatial asymmetry of the potential. For example, it 
seems obvious that the cubic anharmonic Hamiltonian 

d 2 

H = -—r +x 2 +gx 3 (2) 

allows particles to escape to infinity for any real coupling g ^ 0. 

A deeper perturbative and Borel-summation analysis of this problem proved 
such an oversimplified thumb rule wrong many years ago. Calicetti et al j7j 
have shown that the spectrum of the cubic anharmonic oscillator (||) becomes 
purely real and positive at any purely imaginary coupling g — i gi with gi > 0. 
Moreover, the interpretation of the set of resonances of the Hamiltonian (||) as 
analogous to a set of bound states requires a suitable analytic continuation in 
g and/or a careful deformation of the integration path 

The Schrodinger equation for the Hamiltonian operator eq. (Q) cannot be 
integrated algebraically, and one has to resort to numerical methods for its anal- 
ysis. Cannata et al || have found the way of generating complex Hamiltonians 
related to real Hamiltonians, by applying a method invented by Darboux more 
than one hundred years ago. Unfortunately, the resulting models have a very 
complicated structure. 

In this paper we propose a family of problems characterized by two parame- 
ters, see eq. ([|) below, one of them (a) playing a role analogous to the parameter 
a of eq. ([j]) and the other (/3) serving to tune the interaction thus widening the 
richness of the spectrum. Our study of this family will be numerical or, in other 
words, phenomcnological, hoping that it will help a further mathematical anal- 
ysis. One of the main simplifications of our model is that the required complex 
integration paths are lines parallel to the real axis, easy to implement and to 
interpret. 
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2 The model 



We consider the two-parameter family of one-dimensional potentials 

V a p = — (i sinhcc)" cosh' 3 x, (3) 

for arbitrary real values of a and [3. These functions have in general a branch 
point at x — and we select such a branch that the real part of the potential 
is symmetric and the imaginary part antisymmetric with respect to the origin. 
Specifically the potential will be defined by the two equations 

x > V(x) = e m{2+a)/2 sinh Q x cosh x 
x<0 V(x) = e" 27r(2+Q)/2 sinh Q |x| cosh' 9 x 

thus having an invariant Hamiltonian under the PT transformation, i.e., parity 
transformation and complex conjugation. This requires to cut the complex x- 
plane from x = up to x = -co, and to consider the relevant negative x values 
below the cut, i.e. with a small negative imaginary part or a phase — 7r. 

Table 1: Variability of our PT symmetric potential with parameter a, for x > 
and any real j3 



a 






Real part 


Imaginary part 









V R <0 


Vi =0 


< 


a 


< 1 


V R <0 


Vi < 


a = 


1 




v R = o 


Vi < 


1 < 


o 


< 2 


V r >0 


Vi < 


a = 


2 




V r >0 


Vi = 


2 < 


a 


< 3 


V r >0 


Vi > 


a = 


3 




Vr = 


Vi > 


3 < 


a 


< 4 


Vr<0 


Vi > 


a — 


4 




Vr<0 


Vi = 



The characteristic values of the potential for different values of the parameter 
a are shown in Table ||, where we only see the region a 6 [0,4], because the 
same structure is repeated for larger values of a with period 4. In all cases, 
both the real and imaginary parts of the potential tend to either +oo or — oo at 
long distances, when x — > ±oo. The value of f3 does not affect the main features 
of the potential, except for the fact that negative values of /3 may give rise to 
non-confining potentials. We exclude such cases from present study and always 
consider a + (3 > 0. 

The interesting point about PT symmetric Hamiltonians 

H = -f 1 + V aP {x) (4) 



4 



is that they may have localized solutions, corresponding to real eigenvalues, 
which may be interpreted either as bound states or as zero width resonances. 

The special case a — 2 corresponds to a positive and confining potential, 
with an infinite number of bound states. It will be a suitable reference point for 
calculations corresponding to other values of a. We will move towards a > 2 
and a < 2 starting from a given eigenvalue of the case a = 2 to show the 
evolution of the real-energy eigenstates. 

Analogously, one may take any other reference value of a = ckr which pro- 
duces a real and positive confining potential, like or — 6, an — 10, and so 
on. Selecting a different value of an means to consider a different eigenvalue 
problem, which will be labeled by that particular reference value an- Here we 
will concentrate on the case an = 2 mentioned above, and rise considerations 
also for other values of an. 

This family of potentials is similar to the one-parameter family V(x) = 
— (ix) a recently considered by Bender and Boettcher ffl, but the exponential 
growth of our potentials at long distances simplifies the analysis of their prop- 
erties. 

3 Paths in the complex plane 

By carrying out the integration of the Schrodinger equation 

-Z-^ + V^*(x)=E9(x) (5) 

along the real axis (see below for more details), for values of a e [0, 4[, one 
observes that there are solutions for real values of the energy E, smoothly con- 
nected with the solutions of the real potential with a — 2. The wave functions 
are complex, and may be chosen to have a symmetric real part and an antisym- 
metric imaginary part. The last statement is a consequence of the PT invariance, 
which requires that x) is a solution of the Schrodinger equation if ^(x) is 
a solution for a real eigenvalue E. On choosing an appropriate phase factor one 
may then have the mentioned symmetries for the localized wave functions. 

As discussed above, the point a — 2 is an exception, because the Hamiltonian 
is separately parity and time-reversal invariant. Therefore one can choose the 
wave functions to be real, and either symmetric or antisymmetric with respect 
to x = 0. 

At a = 4 the potential is real and everywhere negative; consequently there 
are no eigenstates with real eigenvalues. Real eigenvalues appear again for 
a G]4, 8[, which are smoothly connected with the case a = 6. The same pattern 
is repeated at every a = 2 + AN , for positive integer values of N, as we increase 
a. 

These facts are shown in Fig. |l|, for the case = 0. There, we clearly 
see the existence of two (in general, many) different problems, one centered at 
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Figure 1: Integration of the Schrodingcr equation along the real a;- axis, for the 
model potential V a p. Only the ground-state energy is shown. The two branches 
correspond to solutions smoothly connected with the eigenvalues for aji = 2 
(left curve) and czr = 6 (right curve). The dashed line depicts the expected 
eigenvalues after an appropriate analytic continuation into the complex x-plane, 
for the solutions smoothly connected with = 2. 



E 
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1 .50 .. 




1 .25 .. 



1 .00 



cxr — 2 and the other centered around an — 6. The purpose of this section is 
to find a path in the complex x-plane to obtain lines like the dashed one, which 
represents the continuation of the eigenvalue beyond the limit a = 4. 

The first question to consider is whether there arc confined solutions. To 
this end we take into account the limit x +oo in the Schrodingcr equation, 
where the potential is dominated by the exponential part 

V aP {x) -» exp[*7r(2 + a)/2] exp[(a + (3)x]/2 a+ P . 

One proceeds as in the WKB method assuming a general solution of the form 
^(x) — exp[G(x)]. The leading order of the asymptotic expansion for G{x) is 

where the plus and minus signs come from a square root which appears in the 
differential equation for G(x). 

Except for the particular cases a — 47V, N=0, 1, . . . , there appears to exist a 
solution such that the real part of G(x) is negative and its magnitude increases 
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exponentially, suggesting a discrete set of eigenvalues with localized solutions. 
As mentioned above, we require a + (3 > in order to have asymptotically 
vanishing solutions at long distances. 

The only general statement to be drawn from the above asymptotic limit is 
that, as far as the phase 7r(2 + a)/4 is different from a half integer multiple of 
7r, there are two possible asymptotic solutions, one growing and the other one 
decreasing at long distances; therefore, one one may expect to find one or more 
values of E that select the asymptotically vanishing solutions corresponding to 
localized states. 



3.1 The a R = 2 family 

When a = 2 the phase factor in eq. (||) is exp(i7r), and the required solution 
is the one with the plus sign in front of it (remember that \& = exp[G(x)]). 
For a = 2 + S the phase factor changes to tt + Sir/ 4, and the solution behaves 
asymptotically as 

(a:) exp [[- cos(c57r/4) - i sin(c57r/4)]e (a+/3);r/2 

when x —> oo; i.e., the exponentially decreasing part subsists as far as \5\ < 2. 
The farther we move from a = 2 the slower the exponentially decreasing part 
approaches to zero, indicating that the asymptotic regime will be reached at 
much larger distances. 

There is a way of both moving Asymptotia closer and extending the inte- 
gration beyond a = 4, which consists in adding to the integration variable an 
imaginary component iy which should be negative to be consistent with the po- 
tential branch required by the PT symmetry. After the replacement x <— x + iy 
in eq. (^) the phase angle of the auxiliary function G is changed to 

7r(q + 2) y(a + 0) 

The line corresponding to 9 = tt labels the path of fastest decrease of the 
exponential (the anti-Stokes lines of Ref . || ) , and the boundaries of the region 
of convergence result from the solutions of 9 = 7r±7r/2. 
The optimal path corresponds to a value of y given by 

(2 — o)tt , . 

and remains negative for any value of a > 2, fulfilling the requirement for PT 
symmetry. The integration, however, may be performed for any value of y in 
the range 

(4 - a)ir 

V+ = 2(a + (3) 
air 

V- = 



2(a + /3) 
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Figure 2: The continuous line represents the preferred value of the imaginary 
part y added to the coordinate x, for each value of a and for the family of 
cigcnstates smoothly connected with an = 2. The two dashed lines are the 
boundaries of the acceptable values of y. 



y 




/3=o 



1 2 



a 



2 



so that, for a < 2 the integration should be carried out along the real rc-axis to 
satisfy PT symmetry. The optimal and boundary values of y are plotted in Fig. 
||. We see that without violating PT symmetry, the integration may be carried 
out along the real axis up to a — 4, 

Here our model differs significantly from that considered by Bender and 
Boettcher In their model, eq. (|l|), the integration must be carried out along 
two symmetric sectors, one in the lower-right complex x-plane, and the other 
symmetric with respect to the imaginary axis. The optimal line is given by 
a;exp(z6*), where 9 — ~(a — 2)/(a + 2)(w/2), in such a way that it tends to 
coincide with the negative imaginary axis for large values of a. In our case the 
displacement y remains bounded in the large-a limit, having the value —tt/2. 
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3.2 The families for or = 6 and beyond 

We can carry out an analogous study for the family of potentials connected to 
a = 6, 10, . . . For the case a = 6 one obtains the optimal displacement 

(6 — a)n 
V = 2{a + /3) 

and the lower limit 

(4 - a)ir 
V ~ = 2{a + (3) 
which is suitable for the integration when a > 4. 

To extend the integration to a < 4 we choose an alternative path given by 
the solution with minus sign in eq. (^]), leading to an optimal displacement 
y = -ir(a + 6)/2(a + (3). 

After this discussion we clearly understand the peculiar behavior of the low- 
est real eigenvalue shown in Fig. [I]. The integration having been carried out 
along the real axis, the lowest eigenvalues jump from family to family when 
crossing the points where the potential is purely imaginary. Actually, the way 
of obtaining the eigenvalue by requiring that only the normalizable component 
survives selects the plus or minus sign in eq. ^. In this way, in addition to the 
lack of continuity in the eigenvalues, one may also observe a sudden jump in the 
phase of the wave function (ie, the sign of the imaginary part near + ) when 
crossing each special case mentioned above. 

3.3 The role of the shift y 

It is not difficult to understand the role of the change of variable x — > x + iy if 
we just carry it out explicitly in eq. (|^). The potential-energy function in the 
resulting Schrodinger equation reads 

V cS (x) = V(x + iy), 

so that the new potential with the above mentioned values of y has a dominant 
confining real part and a much smaller imaginary part. 

The actual effect of the transformation is shown in Fig. || for three values 
of a. Particularly impressive is the case a — 4 which originally was a real and 
negative potential, and after the transformation exhibits a dominant real confin- 
ing component. The transformation of the potential guarantees the connection 
through the special points a = AN . 

4 Numerical integration 

The numerical calculation of the eigenvalues of the Schrodinger equation (||) 
with a complex potential and along a complex path, is as simple as in the case 
of real potentials along a real path. 



Figure 3: Real (left) and imaginary (right) parts of the effective potential ob- 
tained by shifting the real variable x to the lower part of the complex plane 
with the optimal y value. The potentials correspond to the phase selection 
appropriate for the family smoothly connected with a = 2. The calculations 
were performed with (3 = for three values of a = 2, 3 and 4, which label the 
corresponding curve 




-4-3-2-10 1 2 3 v 4 -4-3-2-10 1 2 3 v 4 



We have chosen the simplest algorithm which starts by selecting two extreme 
points, X m ; n and X max , at which the wave function is assumed to vanish, and 
discretizing this interval with a uniform integration step h, defined in terms of 
the number of points N as 



N + l 

In order to preserve the PT invariance in the discretization, it is necessary to take 
-Xmin = -^mai- An integer counter k labels the mesh points as xu = X min + kh. 
Approximating the second derivative by the second differences operator, 

dx 2 h 2 

the continuous eigenvalue problem becomes a discrete one given by a symmetric 
and tridiagonal matrix of dimension N and matrix elements 

Hu = + V k 

It is understood that V4 = V(X m i n + kh + iy) in the equation above, where y is 
the appropriate imaginary shift already described earlier. 

The tridiagonal matrix is symmetric, but, contrarily to the case of a real 
potential, it is not hermitean. The roots of the determinant 

D N (E) = det[i? - IE] = 
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give the eigenvalues approximately. The calculation being greatly facilitated by 
the three-point recurrence relation 

D n (E) =D n ^(E)(H m - E) - Hl^D^i, n = 2,3,...,N (10) 

with the starting conditions Dq=1 and D\ = H\\ — E. This recurrence relation 
exhibits the same structure as in the case of a real potential, but not the same 
properties. In the case of a real potential every term of the sequence is real, 
and the set constitutes a Sturm sequence ^ jlQ] . This property allows one to 
devise a simple search algorithm, combining the counting of sign changes of the 
sequence with the bisection method to determine efficiently the eigenvalues. 

This scheme is not appropriate for our case because the potential is complex. 
However, an important feature of our approach is that PT symmetric potentials 
give rise to determinants which are polynomial functions E with real coefficients. 
The proof is quite simple. First one notices that the diagonal matrix elements 
satisfy the rule 

Ha = H N _ i+1N _ i+1 , 

consequence of the PT symmetry. So, the time reversal operation T, which 
corresponds to changing each element of the diagonal by its complex conjugate, 
is equivalent to applying the parity P operation, which corresponds to changing 
every element of the diagonal by its symmetric with respect to the centre of the 
diagonal (this operation is carried out by a similarity transformation with an 
orthogonal matrix having all elements equal to zero with the exception of the 
counter-diagonal elements which are unity). 

In conclusion, the determinant of H — EI is real for real values of E. This 
is true only if X min = — X max , in such a way that only the determinant D N is 
real, but not the terms of the sequence D n with n < N. This feature allows the 
use of the robust bisection method to determine the eigenvalues. 

In almost all cases, we have set the value of a and then determined the 
corresponding eigenvalue. However, in the neighborhoods of the points where 
two real eigenvalues collapse into a pair of complex conjugate eigenvalues, it is 
more convenient to determine the value of a for a given eigenvalue. In any case, 
the method is simple and robust. 



5 The mutual interplay of a, (3 and energies 

Having arrived at the proper way of extending our calculations beyond the 
special points a — AN, our next step is the recomputation of Fig. |l| including 
some excited levels. 

Fig. |] shows results for several levels and for the two sets smoothly connected 
with a — 2 and a — 6. The main feature of this figure is that there is a one-to- 
one correspondence between an eigenvalue with a close to 2 and an eigenvalue 
of the real confining potential with a — 2. The same situation takes place at 
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a = 6, and it is easily proved by means of perturbation theory for a close to 
2,6,... 

Our numerical calculations suggest that there will be real eigenvalues within 
each family with a greater than the reference value aji = 2, and that for smaller 
values of a the real eigenvalues merge into pairs of complex conjugate values 
of E, until reaching the vicinity of a = 1 where once again real eigenvalues are 
allowed. 

Figure 4: Spectrum of the PT invariant Hamiltonian with the potential V a p 
and ft = 0, showing several bound states corresponding to two families, one 
connected with a = 2 (continuous line) and the other connected with a = 6 
(dashed lines) 



2 

20 




I | | | | | | | | | 

1 2 3 4 5 6 7 8 9 10 

Figure || shows in detail the special characteristics of the levels near a = 1 for 
ft = 0. In particular, this figure illustrates the simultaneous jump of the fourth 
and fifth levels into the complex plane, near a = 1.15, and their simultaneous 
return to the real axis when a is slightly greater than unity. The same pattern 
seems to happen also for higher levels. Obviously, as far as the characteristic 
polynomial is real, the transition from real to complex values must be in pairs. 
We have not observed a similar phenomenon in the case of the set connected 
with a = 6, but it may well happen for levels of energy higher than those shown 
in Fig. |. 

Up to now we have concentrated on calculations with = 0. Fig. || illus- 
trates the role of the parameter j3. In addition to the case j3 = 0, this figure 
displays also the lowest levels for several values of ft, both positive and negative. 
The interesting role of ft is to switch the special level, i.e., the level which ulti- 
mately will move around a = 1. With its help one may choose this special level 
to be the first one, (ft = -0.25), the third one (ft = 0), the fifth one (ft = 0.25) 
and so on. Because of the special way the levels form the pairs, one should 
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Figure 5: Enlarged view of the spectrum for j3 = near a = 1. 
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not be surprised that even levels cannot become special in the above mentioned 
sense. 

6 Quasi-algebraic study 

As in our previous paper on the subject, we have also supplemented the nu- 
merical calculation with the Riccati-Pade method (RPM). When a = 2 and 
(3 = the potential-energy function is parity-invariant and the RPM leads to 
just one Hankel determinant from which one obtains the eigenvalues jl2|]. The 
calculation is straightforward and the rate of convergence sufficiently great as 
shown in Table ||. 

Table 2: The RPM ground-state energy for j3 = and a = 2 in terms of the 
dimension of the Pade determinant 



hline D 


RPM root 


2 


1.213616523 


3 


1.211409311 


4 


1.211411109 


■5 


1.2114109830 


G 


1.211410984169 


7 


1.2114109841755 


8 


1.21141098417527 
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Figure 6: Several eigenvalues computed at f3 = 0.5, (3 = 0.25, (3 = and 
= -0.25 
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For nonhcrmitean cases we change the coordinate according to x — iq so 
that the Hamiltonian operator becomes 

and we can apply the RPM as in the case of a real Schrodinger equation. If the 
potential energy is an even function of q we apply the method just indicated; 
if it is not then the RPM leads to two Hankel determinants E3] from which we 
obtain both E and -i*'(0)/*(0). Table | shows results for a = 1 and a = 3 
in excellent agreement with the numerical integration discussed above. 

We have carried out the RPM calculations algebraically by means of Maple, 
resorting to a numerical approach just at the end in order to obtain the roots of 
the Hankel determinants |L2| [l3| . For this reason the requirement of computer 
memory is considerable in the case of a nonsymmetric potential-energy function, 
and we cannot handle determinants of the same dimension as in the symmetric 
case. For (3 = we have tried to overcome this problem by means of the change 
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Table 3: The RPM ground-state energies for the non-hermitean cases a = 1 
and a = 3. In both cases is (3 = 0. 





a = 1 


a = 1 


a = 3 


a = 3 


D 


E 


z*'(0)/*(0) 


E 


i*'(0)/*(0) 


2 


1.655005966 


1.033573034 






3 


1.765033153 


1.095023981 


1.385656774 


-0.5049697062 


4 


1.765157398 


1.095137449 


1.349869536 


-0.4769952880 


5 


1.765157246 


1.095137384 


1.350149473 


-0.4771536171 


6 






1.350140759 


-0.4771520606 


Numeric 


1.76515725 


1.09513737 


1.350140990 


-0.47715200 



Table 4: The RPM ground-state energies for the non-hermitean cases a = 1 
and a = 3 with the complex rescaled Hamiltonian. In both cases is = 0. 





£(a = 1) 


E(a = 3) 


2 


1.765248635 


2.60904086409 


3 


1.765157328 


2.59510727493 


4 


1.765157255 


2.59524841637 


5 


1.76515725525231 


2.59524599823 


6 


1.76515725525336 


2.59524605087 


7 


1.7651572552533587 


2.59524605034 


8 


1.76515725525335874 





of coordinate x = i(u + n/2) jnj and applying RPM for symmetric potential- 
energy functions to the resulting Hamiltonian operator 

d 2 

-H = -— + [-cos u] a . 

Table || shows results for a = 1 and a = 3. In the former case the result is 
identical (though more accurate) to the one shown in Table ||; however, in the 
latter case we do not obtain the lowest eigenvalue but the first excited one. We 
have not yet being able to justify this anomalous behavior of the RPM when 
a = 3. 
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